pmf

We can use the base plot functions in R to create a plot of the pmf for a binomial random variable \(X\) with \(n=13\) and \(p=0.6\) — i.e. \(X\ \tilde \ B(13,\ 0.6)\).

  n <- 13
  p <- 0.6
  x <- -1:14
  pmf <- dbinom(x, n, p) ### dbinom() is the pmf
  rbind(x, pmf)
##     [,1]         [,2]         [,3]        [,4]        [,5]       [,6]
## x     -1 0.000000e+00 1.0000000000 2.000000000 3.000000000 4.00000000
## pmf    0 6.710886e-06 0.0001308623 0.001177761 0.006477683 0.02429131
##           [,7]      [,8]      [,9]     [,10]     [,11]      [,12]       [,13]
## x   5.00000000 6.0000000 7.0000000 8.0000000 9.0000000 10.0000000 11.00000000
## pmf 0.06558654 0.1311731 0.1967596 0.2213546 0.1844621  0.1106773  0.04527707
##           [,14]        [,15] [,16]
## x   12.00000000 13.000000000    14
## pmf  0.01131927  0.001306069     0
  cbind(x, round(pmf,5))
##        x        
##  [1,] -1 0.00000
##  [2,]  0 0.00001
##  [3,]  1 0.00013
##  [4,]  2 0.00118
##  [5,]  3 0.00648
##  [6,]  4 0.02429
##  [7,]  5 0.06559
##  [8,]  6 0.13117
##  [9,]  7 0.19676
## [10,]  8 0.22135
## [11,]  9 0.18446
## [12,] 10 0.11068
## [13,] 11 0.04528
## [14,] 12 0.01132
## [15,] 13 0.00131
## [16,] 14 0.00000
  plot(x, pmf, type="h", xlab="x", ylab="p = P(X=x)", main="X~B(13, 0.6)", 
       ylim=c(0, 0.25), xaxt="n")
  axis(1,at=seq(0,14,by=2))
  text(x, pmf+0.005, round(pmf, digits=4), cex=0.6)

We note that the maximum is at 8 .

CDF

The CDF may be plotted analogously.

  x <- -1:15
  cdf <- pbinom(x, n, p) ### pbinom() is the CDF
  plot(x, cdf, type="s", xlab="x", ylab="P(X<=c)", 
       main="X~B(13, 0.6)", xaxt="n")
  axis(1, at=seq(0,14,by=2))

Just for fun we can overlay the two.

  pmf <- dbinom(x, n, p)
  plot(x, pmf, type="h", xlab="x", ylab="p", main="X~B(13, 0.6)", 
       ylim=c(0, 1), xaxt="n", col="red", lty=3)
    lines(x, cdf, type="s", xlab="x", ylab="P(X<=c)", 
          main="X~B(13, 0.6)", xaxt="n")
  axis(1,at=seq(0,14,by=2))

Most Probable Value

Looking at successive terms allows us to find the most probable values of \(X\). We note that \[ \frac{p(x+1)}{p(x)} = \frac{n-x}{x+1}\cdot\frac{p}{1-p} \] We can look at where the ratio is greater/less than one.

  x <- 0:n
  succ_terms = ((n-x)/(x+1))*(p/(1-p))
  plot(x, succ_terms)
  abline(h=1, lty=2, col="green")

  cbind(x, succ_terms)
##        x succ_terms
##  [1,]  0 19.5000000
##  [2,]  1  9.0000000
##  [3,]  2  5.5000000
##  [4,]  3  3.7500000
##  [5,]  4  2.7000000
##  [6,]  5  2.0000000
##  [7,]  6  1.5000000
##  [8,]  7  1.1250000
##  [9,]  8  0.8333333
## [10,]  9  0.6000000
## [11,] 10  0.4090909
## [12,] 11  0.2500000
## [13,] 12  0.1153846
## [14,] 13  0.0000000

Another way of looking at the problem is to compute the ratios by using the pmf.

  pmf <- dbinom(x,n,p)
  cbind(x, pmf)
##        x          pmf
##  [1,]  0 6.710886e-06
##  [2,]  1 1.308623e-04
##  [3,]  2 1.177761e-03
##  [4,]  3 6.477683e-03
##  [5,]  4 2.429131e-02
##  [6,]  5 6.558654e-02
##  [7,]  6 1.311731e-01
##  [8,]  7 1.967596e-01
##  [9,]  8 2.213546e-01
## [10,]  9 1.844621e-01
## [11,] 10 1.106773e-01
## [12,] 11 4.527707e-02
## [13,] 12 1.131927e-02
## [14,] 13 1.306069e-03
  succ_pmfs <- pmf[2:(n+1)]/pmf[1:n]
  succ_pmfs <- c(succ_pmfs,0) ### Last term is 0 because p(n+1) = 0
  cbind(x,succ_pmfs)
##        x  succ_pmfs
##  [1,]  0 19.5000000
##  [2,]  1  9.0000000
##  [3,]  2  5.5000000
##  [4,]  3  3.7500000
##  [5,]  4  2.7000000
##  [6,]  5  2.0000000
##  [7,]  6  1.5000000
##  [8,]  7  1.1250000
##  [9,]  8  0.8333333
## [10,]  9  0.6000000
## [11,] 10  0.4090909
## [12,] 11  0.2500000
## [13,] 12  0.1153846
## [14,] 13  0.0000000

Plotting the ratios generates the same plot as above.

  x <- 0:n
  plot(x, succ_pmfs)
  abline(h=1, lty=2, col="green")

In either case, we see that the function goes from increasing to decreasing at the integer part of \((n+1)p\). This can be computed as \(\mbox{trunc}((n+1)p)=\) 8 or \(\lceil np \rceil =\) 8 .

  (max_x <- trunc((n+1)*p))
## [1] 8
  (max_x <- ceiling(n*p))
## [1] 8